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LONG-TIME NUMERICAL INTEGRATION OF THE THREE-DIMENSIONAL WAVE 
EQUATION IN THE VICINITY OF A MOVING SOURCE* 

V. S. RYABEN’KlU, V. I. TURCHANINOV* , AND S. V. TSYNKOV§ 


Abstract. We propose a family of algorithms for solving numerically a Cauchy problem for the three- 
dimensional wave equation. The sources that drive the equation (i.e., the right-hand side) are compactly 
supported in space for any given time; they, however, may actually move in space with a subsonic speed. 
The solution is calculated inside a finite domain (e.g., sphere) that also moves with a subsonic speed and 
always contains the support of the right-hand side. 

The algorithms employ a standard consistent and stable explicit finite-difference scheme for the wave 
equation. They allow one to calculate the solution for arbitrarily long time intervals without error accumu- 
lation and with the fixed non-growing amount of the CPU time and memory required for advancing one time 
step. The algorithms are inherently three-dimensional; they rely on the presence of lacunae in the solutions 
of the wave equation in oddly dimensional spaces. 

The methodology presented in the paper is, in fact, a building block for constructing the nonlocal highly 
accurate unsteady artificial boundary conditions to be used for the numerical simulation of waves propagating 
with finite speed over unbounded domains. 

Key words, wave equation, lacunae, finite-difference approximation, explicit numerical integration, 
arbitrarily long time intervals, non-accumulation of error, uniform error bounds, fixed expenses per time 
step 

Subject classification. Applied and Numerical Mathematics 

!• Introduction. We will be solving numerically a Cauchy problem for the three-dimensional wave 
equation: 


d 2 u 
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/ d 2 u d 2 u & 2 u\ 
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We will bp interested in calculating the solution u = u(x,t) only for those values of the argument x = 
(j’i . x->, x :() that belong to the ball S = S(t) defined by the following inequality 


(xi - x?) 2 + (x 2 - x? 2 )' 2 + (x 3 - x!>) 2 < d 2 / 4. (1.3) 

In inequality (1.3), d is the diameter of the sphere S(t), and the functions x° = x^t), x!J = x° 2 (t), and 
^3 = •'■•((O are given smooth functions of their argument t — time; these functions define the motion of the 
sphere’s center so that 



(1.4) 


where k < c. In other words, we consider only subsonic motions of the sphere S(t ), the maximum speed 
k of its center may never exceed the speed of sound c (see equation (1.1)). Regarding the right-hand side 
f(x,t) of equation (1.1), we always assume that it is a sufficiently smooth function with respect to all its 
four arguments and also 


suppf(x,t) = {(xi,X2,ar 3 ,0l® € S{t), t > 0} . (1.5) 

In this paper, we construct efficient numerical algorithms for solving the foregoing problem. These 
algorithms employ the simplest central finite-difference scheme that approximates problem (1.1), (1.2) on 
smooth solutions with the second order of accuracy. This scheme has been chosen primarily for the reason 
of simplicity so that we can demonstrate the concept using the least cumbersome approach. Generally, any 
scheme that possesses the properties of stability and consistency on smooth solutions, including high-order 
schemes, can be used for building the algorithms of the type described hereafter. 

The second order central-difference scheme is constructed on the uniform Cartesian grid with the size h 
in all spatial directions and time step r: 

j(^mi z ma jtm 4 ) = (m t ft, m 2 /i, m 3 ft, m 4 r) |m i , m 2 , m 3 = 0, ±1, ±2, . . . , rn 4 = 0, 1,2, ... J. 

In every grid node ni = (raj , m2, m3, 7714), equation (1.1) is replaced by the finite-difference equation 


^ ^ Q'mri'Mn — frm — 2 , 3, 4 , . . . . ( 1 . 6 ) 

The discrete solution u n = u nitWa , n3i n 4 is defined on the same grid: ni,n 2 ,n3 = 0, ±1,±2, ..., n 4 = 
0, 1,2,3,...; equation (1.6) connects the values of u n in the following nine nodes of the grid that form 
the stencil N m : 


(mi/i,m 2 /i, m 3 /i, (m 4 - j)r), j -0,1,2, 
((mi ± l)/i, m 2 /i, m$h, (m 4 - 1 )r), 

(5 ni\h , (m 2 ± 2 )h,m 3 h, (m 4 - 1 )r), 
(mi/i,m 2 /i,(m3 ± 1 )h,(m 4 - l)r). 
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The coefficients a mn of the discrete operator and the values f m of the discrete right-hand side in (1.6) are 
defined as follows 


1, 
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fm = (m 4 - 1 )rh), (1.9) 

where r is the Courant number 


r = 


r 

- < 
h ~ 


1 

cy/i' 


( 1 . 10 ) 


estimate (1.10) immediately follows from the standard stability considerations of von Neumann type. Initial 
conditions (1.2) are replaced by the conditions 


«?.,,n a ,n 3 = °> p = n 4 — 0, 1. (l.H) 

It is known that if the right-hand side f(x,t.) is sufficiently smooth with respect to all its arguments , then the 
finite- difference solution 


U n — ?im ,ni>,n 3 ,n 4 — w nl.ti2,«3’ 

where p — n 4 , converges to the continuous solution u(x,t) of problem (1.1), (1.2) as the grid size decreases , 
and the following estimate holds 


max \u(nih,n 2 h,n 3 h,pT ) - , n2>n3 | < C\f\ KiT h'\ pr < T fina |, (1.12) 

n 1 ,»2 

here C = <7(7^1) is a constant, A' is a sufficiently large positive integer number, and \f\x,T I s the sum of 
maximal absolute values of all derivatives of the function f(x,t) up to the order K for t < TH na i. Therefore, 
speaking formally, estimate (1.12) allows one to use the finite-difference scheme (1.6), (1.11) for calculating 
the solution u(x, t) of problem (1.1), (1.2) on arbitrarily long time intervals 0 < t < Tfinai. 

There are, however, two most substantial obstacles. First, when calculating the solution using the scheme 
(1.6), (1.11), the number of grid nodes involved in the computation on each time level increases approximately 
as (d/h + p ) 3 with the number of level p. Consequently, when p & Thnai/r and r = r/i, r - l/c\/3, this 
number is of the same order of magnitude as Therefore, the associated storage and CPU time 

requirements grow rapidly as the final time T fina i increases. Second, estimate (1.12) basically does not 
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guarantee the algorithm from the possible error accumulation as the final time T final and the corresponding 
number of steps increase. Indeed, although the constant C in inequality (1.12) does not depend on the grid 
size ft, it, generally speaking, depends on the final time 7fi na i (may grow with T fina |). 

In the paper, we propose three approaches to improve the simplest scheme (1.6), (1.11); these approaches 
take into account some specific properties of the solution to problem (1.1), (1.2). Similar approaches can 
be developed for other converging schemes, besides the scheme (1.6), (1.11). Each of the three algorithms 
proposed hereafter guarantees that the error will not accumulate as the number p of the time level in- 
creases. Moreover, both the memory and CPU time required for advancing each time step remain bounded 
independently of p (and Tk na j) for the fixed grid sizes ft and r. 

The number of arithmetic operations required for advancing one time step using either of the three 
algorithms that we are discussing is 0(A r ), where N is the number of grid nodes in space (i.e., on one time 
level) inside the sphere of the fixed diameter d\ clearly, N — 0(h~ 3 ). This number does not depend on t, i.e., 
does not increase with p because unlike in the original scheme (1.6), (1.11), the computational domain in the 
new algorithm will not need to expand in space as the time goes by. Obviously, the number O(N) is optimal 
(linear with respect to the grid dimension) and cannot be improved by choosing any other algorithm. The 
required memory (number of words) in all three versions of the new algorithm is of the order 0{N) as well, 
i.e., does not exceed some C N, where C is a constant that depends on the version. For the third algorithm 
this constant is smaller than for the others. 

All three aforementioned methodologies for improving the original scheme (1.6), (1.11) so that one can 
calculate the solution u(x,t) of problem (1.1), (1.2) on the domain x € S(t) and arbitrarily large time 
intervals rely on a particular property of solutions to the three-dimensional wave equation (1.1), namely the 
property of having lacunae. 

The lacunae-based technique proposed in this paper for the long-time numerical integration of the three- 
dimensional wave equation in the neighborhood of a moving source is a generalization and extension of the 
technique developed previously in [1] for the case of stationary sources. In the future, these lacunae-based 
techniques will be used as a building block for constructing global artificial boundary conditions (ABCs) for 
the numerical simulation of waves propagating with finite speed over infinite domains. The latter framework 
includes, in particular, the problems of electromagnetic diffraction and scattering, as well as the problems in 
both ambient and advective acoustics. The unsteady ABCs’ methodology that we have mentioned will be 
described in a forthcoming publication. A general review of different ABCs’ methodologies available in the 
literature can be found in [2]. 

2. Lacunae of the Wave Equation. Consider the nonhomogeneous wave equation with zero initial 

conditions: 


d 2 v 2 / d 2 v d 2 v d 2 v \ 
dt 2 C \ dx{ + dx 2 dx 2 ) 


t > 0 , 


(2.1) 


Assume that <p(x,t) is a sufficiently smooth function with respect to all its arguments and that ip{x,t) = 0 
for t < 0. 
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For every (;r, t), the solution v = v(x, t) of problem (2.1), (2.2) can be written in the form of the Kirchhoff 
integral: 


v(x,t) 



p<ct 


p 


d£, 


(2.3) 


where p = \/(x\ — £i ) 2 + ( x 2 — £j) 2 + (#3 — £3)-' and d£ = d^id&dfa. The integration in (2.3) is performed 
over the ball of radius ct centered at x in the space £ = (£1 ,£>,£;i)- Formula (2.3), in fact, implies that the 
solution v{x , t) at the point (x, t ) depends only on the values of ip(£, 6) on the surface of the characteristic 
cone (its lower portion) with the vertex 


(zi - 6) 2 + (*2 -6) 2 + (^3 -6) 2 = c\t-0f, 6 <t (2.4a) 

and does not depend on y>(£, 9) when (£, 9) belongs to the interior of the cone (2.4a). In other words, changing 
the values of p(£,9) in the interior of the cone (2.4a) will not affect the solution v{x, t) at the point ( x,t ). 
To emphasize this circumstance, we will call the domain 


(^i _ £i) 2 + (^2 — $i) 2 + (%3 ~ <^3) 2 < c 2 (t - 0 ) 2 , 9 < t, (2.4b) 

i.e., the interior of the characteristic cone (2.4a) the lacuna of the right-hand side of equation (2.1) with 
respect to the point (x,t,). The presence of the lacuna (2.4b) of the right-hand side implies that the solution 
v(x,t) of (2.1), (2.2) will also have a lacuna D+(Q) with respect to the domain Q of the right-hand side. 
Indeed, consider a 6-source for equation (2.1) concentrated at the point (f,0) of the space-time: s&e). At 
any moment of time t > 9, the solution of problem (2.1), (2.2) driven by this source will be concentrated on 
the surface of the sphere of radius c(f - 9) centered at £ in the space x = (xi , x 2 , x 3 ). Inside this sphere, the 
solution will be identically zero: v(x,t) = 0 for p = \x — £| < c(t — 9). Therefore, let us now interpret the 
surface (2.4a) as the upper portion of the characteristic cone of equation (2.1) in the space-time (x, t) with 
the vertex (£, 9). Then, the solution of (2.1), (2.2) driven by 6(£, 9) is zero in the interior of the cone (2.4a), 
i.e., on the domain (2.4b) that we now denote by D + (£, 9), 


D + {l0) = {{x,t)\\x-t\<c{t-8), t>8 }, (2.4c) 

and call the lacuna of the fundamental solution of the wave equation. (Note, this fundamental solution is 
actually a single layer on the spherical surface |x-£| = c(t-9 ), t > 9.) If we consider a general source ip(£, 9) 
rather than the 6-source 6(£, 0), then for every particular (£, 9) the solution of (2.1), (2.2) inside the lacuna 
£>+(£, 0) given by (2.4c) does not depend on the value of y>(£,0) at this point (£,0). By the superposition 
principle, the solution of (2.1), (2.2) with a general source ip will be concentrated on the set given by the 
union of all spheres \x - £| = c(t - 0), t > 0, when the vertex (£,0) of the cone (2.4a) sweeps the support of 
the right hand side supp <p(£,9). Accordingly, the intersection of all £>+(£, 9) of (2.4c) for all (£, 0) e Q will 
be called the lacuna of the solution v(x,t) with respect to the domain Q\ 
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(2.4d) 


D + (Q)= fj D + {1,6). 

U,B)€Q 

Clearly, the solution v(x,t) of (2.1), (2.2) is zero on D + (Q) of (2.4d), 

v(x,t) for (x,t) € D + {Q) (2.5) 

if 


supp^C Q. (2.6) 

Alternatively, one can say that changing the values of 6) in the domain Q is not going to affect the 
solution n(x,f) of (2.1), (2.2) in the points of the lacuna D + (Q) given by (2.4d). In other words, we see 
that the waves governed by the three-dimensional wave equation (2.1) have trailing fronts. If the source 
is compactly supported in both space and time, then at any given location x in space the solution v(x.t) 
becomes identically zero after a finite interval of time. This finite time interval is the time from the moment 
of source activation till the moment when the point x falls into the lacuna D + (Q) given by (2.4d), or in 
other words, till all the waves generated by the source have pass through x and accordingly, the solution 
there become zero again. 

If the domain Q is defined as follows 


Q = {(x,t)|x € S(f), t 0 <t<t a, 


(2.7) 


then condition (2.G) implies that that the solution v(x,t) of (2.1), (2.2) satisfies the identities 


and 


v(x, t) = 0, for t < to 


(2.8a) 


v(x,t) = 0, for X € 5(f), t > t 0 + - - - _*° fc )(C + . (2.8b) 

The first identity, (2.8a), is obvious, it takes place because the initial data of the Cauchy problem are 
homogeneous, see (2.2). The second identity, (2.8b) holds in virtue of (2.6) because the region of the space- 
time ( x,t ) defined as (x € S (t), t>t 0 + d +( f| zis!L£±tl j j s completely contained inside the lacuna D + (Q ) 

of (2.4d). In other words, as long as (1.4) holds the time interval , s sufficient for all the waves 

generated by the sources inside S(t ) during t 0 < t < t j to completely leave the moving domain S(t). For 
the case of stationary sources, k = 0, the inequality of (2.8b) reduces to the obvious estimate t > tj + d/c, 
set' [1], Let us also note that the estimate for t given in (2.8b) is, in fact, conservative, it does not make 
any assumptions regarding the character of the source motion except that its maximal speed is k < c. To 
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obtain (2.8b), one only assumes that at any t > to the source, i.e., 5(f), can be anywhere inside the sphere 
of diameter d + 2 k(t — fo) centered at the center of 5(fo). If, however, we make an additional assumption 
regarding the motion of the sphere, e.g., that it moves with a constant speed k is some prescribed direction, 
then the estimate of (2.8b) can be alleviated and instead we obtain 


w(a;,t)=0, for x e 5(f), f > t\ H 

c — k 

For a stationary source, k — 0, (2.8c) again reduces to f > f] + d/r of [1]. 

Let us now introduce the following partition of unity. Define the function 


(2.8c) 


- < 


1, 

0 < t < 1/2 

m 

1/2 < t < 1 

o, 

t > 1 

*(-o, 

t < 0 


where P(t) is a, polynomial of the type 


(2.9) 


with the coefficients a, b, c, d, and e such that the following equalities hold 


i J (t) = I + a (*-|')+6(*-7 


P( 1) = P # (l) = P"(l) = P'"( 1) = P iIV) ( 1) = 0. 

Obviously, P(l/2) = 1 and the derivatives of P(t) up to the fourth order are equal to zero at t = 1/2. 
Therefore, 'f'(f) is an even function with four continuous derivatives for all t € R and also 'F(f) is compactly 
supported, ^(f) = 0 for |f| > 1, i.e., 


supp^(f) = [-1, 1]. 

Specify now a parameter T and introduce the functions 


9j{t,T) = 9 , j = 0,1,2,... 

Clearly, 


supp = [O' - 1 )T, 0 + 1 )T], j = 0,1,2,-.. 


Moreover, for any T > 0 we have 
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( 2 . 10 ) 


oc 

'52*j(t,T) = l, t> 0. 

j = o 

The representation of a function which is identically equal to 1 in the form (2.10) is a partition of unity. 
Notice that for every given t no more than two terms on the left-hand side of the identity (2.10) may differ 
from zero. 

We now represent the right-hand side f(x,t) of equation (1.1) in the form 

oo oc oc 

f(x,t) = f(x,t)^j(t,T) = T)f(x,t) = ^Tfj(x,t,T), (2.11) 

J=0 j= 0 j= 0 

where fj(x,t,T) — j(t,T)f(x,t). Clearly, 


supp fj{x,t) = Qj(T), 


( 2 . 12 ) 


QAT) = { (*,t)l * e S(t), (j-i)T<t< u + i)rj . 

Consider the following sequence of problems 


(2.13) 


d ~ u i _ 2 ( q2u j , d 2 uj 

dx 2 Ox 2 dx 2 


di 2 


j = fj{x,t,T), 


_ dUj_ 
#=o dt 


(2.14) 


= 0, j = 0,1,2,... 


/ = () 


Because of the linearity of problem (1.1), (1.2) and representation of f(x,t) in the form (2.11), the solution 
u(x,t) of problem (1.1), (1.2) can also be represented as a similar sum 


oc 

u{x,t) = £uj(x,t,T), (2.15) 

j=o 

where uj(x,t,T) is the solution of problem (2.14). Let us show that for x £ S(t) and any fixed t > 0 there 
are only a few values of j for which uj(x,t,T) / 0. First, we apply identities (2.8a) and (2.8b) which hold 
under conditions (2.6), (2.7) to the solution uj(x,t,T) of problem (2.14). In so doing, instead of (2.6), (2.7) 
we use (2.12), (2.13). Then, instead of (2.8a) and (2.8b) we obtain the following two identities 


and 


Uj(x,t,T)= 0, for t < (j — \)T 


(2.16a) 


uj(x,i,T)= 0, for x 6 5(f). t > (j-l)T+ — - 

c — k 


(2.16b) 
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Identities (2.16) imply that for any given t and T the solution Uj(x,t,,T) may differ from zero for x € S(t) 
only if the following two inequalities hold simultaneously 


t<(j-l)T + 


U ~ 1 )T < t, 

d + 2 T{c + k) 


(2.17a) 


(2.17b) 


A fixed prescribed t = t° can meet both conditions (2.17) if and only if the index j satisfies the inequalities 


t° d + 2T(c+k) 
T (c - k)T 


<J < 1 + y. 


(2.18) 


Therefore, there is only a finite number of values j for which Uj(x,t,T) differs from zero for x 6 S{t) and 
given t. If k — 0 and T > d/e (and also t° is sufficiently large) then there is no more than three such values 
of j. If T — > +0 or k — > c then the number of indexes j that satisfy (2.18) increases with no bound. 

Henceforth, we will be using representation (2.15) for the solution u(x, t) of problem (1.1), (1.2). We 
note that the term Uj(x,t,T) in formula (2.15) is of interest for us only till the moment 


t = (j-l)T + 


d + 2T(c+k) 
c — k 


(2.19) 


as starting from this moment the component Uj(x,t, T) turns into zero inside the computational domain 
S(t.) because of (2.16b) and therefore no longer contributes into the sum (2.15). The sphere S(t) of diameter 
d centered at (xj {t), x^it), ^(J)) represents at the time t of (2.19) the trailing front of the propagation of 
Uj(x,t,T) over the unperturbed zero background. (In fact, in many cases the spherical surface S(t.) may be 
a conservative estimate for the actual location of the trailing front; but at any rate, S(t) is always inside the 
trailing front.) 

Numerical algorithms proposed hereafter are based on the concept that when calculating the solution 
u(x,t) of (l.l)j for every t we actually need to calculate only a few terms nj(x,t,T) in the sum (2.15) 

that differ from zero for x £ S{t). 

Let us make the following important remark. Specify some z > 0 and consider the following problem 
periodic with the period z in every coordinate direction xi, l = 1, 2, 3. 


d'S _ -i ( , &vi , \ 

dt. 2 \ dx\ dx 2 dx\ J 


fj(x,t,T,z), 


Vj (x, t, T,z) = 0, t < (j - 1 )T, 

Vj(x 1 + kiZ,X-2 + k-2Z,X 3 + k3Z,t,T,z) = t >j(x,t,T,z) 
fj(x i + k[Z,x-2 + k 2 z,x 3 + k^z,t,T,z) = fj(x,t,T,z) 


fci , k'i = 0, ±1, i2, . . . 

fj(x,t,T,z) = fj(x,t,T) if |x(| < 2/2, / = 1,2,3. 


(2.20a) 


(2.20b) 
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THEOREM 2.1. The solution vj(x,t<T, z) of problem (2.20) coincides with the solution Uj(x,t,T) of 
problem (2.14) 


Uj(x,t,T) = Vj(x,t,T,z) 

on the domain 


( 2 . 21 ) 


xeS(t), (j-l)T<t<(j-l)T + ^-. 


( 2 . 22 ) 


Proof Let us first note that as long as 


dx°(t) 

dt 


< k < 1 (see (1.4)), where a* 0 = x°(£) may be any 


prescribed law of motion for the center of the sphere S(t ), the right-hand side /j(x, t,T y z) of (2.20a), which 
is periodic in all three coordinate directions xi, x 2 , and x 3 , may differ from zero only on the collection of 
balls 


(Xi 


k[z) 2 + (x 2 - h>z) 2 + (x 3 - k 3 z) 2 < 




t>{j~ 1)7\ *i, k 2 , k 3 =0,±1,±2, ... 


(2.23) 


This actually follows from the fact that the sphere S(t) for t > (j - 1)T completely belongs to the ball (2.23) 
for A - 1 = k 2 = k 3 = 0. Moreover, it is easy to see that the lower portion of the characteristic cone (2.4a): 


(*i - 6) 2 + (*2 - &) 2 + (x 3 - 6) 2 = c 2 (t - 6 < t, 

with the vertex in an arbitrary point (x, t) that belongs to the ball (2.23) for A’, = k 2 = fc 3 = 0 intersects 
none of the spherical domains (2.23) for A-jf + k' 1 , + k% / 0 (i.e., none of the other balls (2.23)) on the time 
interval (j — 1)7’ < 0 < f if only 


t<(j-l)T+ (2.24) 

Then, the Kirchhoff formula (2.3) implies that the value of the solution Vj(x,t,T,z) in the vertex (x, t) of 
the characteristic cone (2.4a) will not depend on the values of the right-hand side fj(£,0,T, z) of equation 
(2.20a) on the domains (2.23) for kf + + & 3 ^ 0 (£ is substituted for x and 0 is substituted for t in formula 

(2.23) ). In particular, the value vj(x, t , T, z) will not change if the right-hand side fj(£, 0, T, z) on all domains 

(2.23) except the “central” one k\ — h 2 = k 3 — 0 was replace by the identical zero for all 6 < t, t of (2.24). 
On the other hand, this replacement makes the right-hand side of (2.20a) coincide with the non-periodic 
right-hand side of equation (2.15), which has the solution Uj(x,t,T). Thus, Vj(x,t,T, z) — uj(x,t,T) for all 
those (x, t) y for which t satisfies (2.24) and x belongs to the ball (2.23) for k\ = k 2 = fc 3 = 0. At the same 
time, it has been mentioned that the sphere S(t) completely belongs to the ball (2.23), k x = k 2 = Jfc 3 = 0, 
for any t > (j — l)T. This proves the theorem. □ 
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3. Finite-Difference Algorithms. We will construct three algorithms for the approximate calculation 
of the solution to problem (1.1), (1-2) on arbitrarily long time intervals using finite-difference scheme (1.6). 
All three algorithms will guarantee that the error will not accumulate and computer expenses per time step 
(both CPU time- wise and storage- wise) will not increase, i.c., will remain fixed and bounded throughout the 
entire computation. All three algorithms will also have the same non-irnprovable computational complexity, 
i.e., asymptotic order of the number of required arithmetic operations and amount of memory with respect 
to the grid size. However, the algorithms will differ from one another by the actual number of required 
arithmetic operations and amount of memory (while still having the same asymptotic order with respect to 
the grid size /i), as well as by certain convenience features from the standpoints of their justification on one 
hand and practical computing on the other hand. 

For all three algorithms we assume that T > 0, h > 0, r > 0, and z > 0 are chosen so that T = qr , 
z — rh, where q and r are positive integers, and r/h < l/c\/3, i.e., the von Neumann stability criterion is 
met. 

The First Algorithm . This algorithm is based on the representation (2.15) of the solution u(x,t) to 
problem (1.1), (1.2) in the sphere S(t). Let us fix some arbitrary integer / and consider t from the interval 


{l-l)T <t<lT. (3.1) 

For these t , formula (2.15) can be rewritten as follows 


u(x,t) = ui- s (x,t,T) + u/_ s+1 (:r,UT) + . . . + ui(x,t,,T), (3.2) 

where Uj(x , T), j = l — / — s + 1, . . . , /, are solutions of the corresponding problems (2.14). The positive 

integer number s is chosen from the inequalities (3.1) and (2.18). If, instead of the smallest possible s that 
satisfies the foregoing constraints, one takes, e.g., s + 1, then an additional term (x, t, T) will simply 

appear in the sum (3.2). This term, however, will turn into zero for x £ S(t ) and t satisfying (3.1) and 
consequently, the work required for computing this term will be superfluous. 

Assume, for definiteness, that s is chosen according to the formula: 


8 = 


d + 2 T(c + k) 
(c-k)T 


(3.3) 


where [ ■ ] denotes the integer part. We will be calculating the solution u(x,t) of problem (1.1), (1.2) on the 
grid n = (n\h, n 2 /i, n 4 r) for 


( n\h , n 2 h, n^h) £ S(n 4 r), (/ - 1)T < n 4 r < IT , 

i.e., inside the sphere S(t) for t satisfying (3.1). According to (3.2) and (3.3), the exact values of this solution 
on the grid are given by 


u(nih , n 2 h, n 3 /i, n 4 r) = u/_ s (t?i/z, n 2 /i, n 3 h, n 4 r, T) + 


-hu/_ s+ i (rii h, n 2 h, n 3 /i, n 4 r, T) -F . . . -b u/(n i h , n 2 h , n 3 /i, n 4 r, T) 


(3.4) 


li 



Instead of the exact values (3.4) the first algorithm that we are discussing yields approximate values of the 
solution u(x, t) on the grid according to the formula: 


u(n 1 h, n- 2 h , n 3 /t, n 4 r) « u\P = u n (l - s,T) + u n {l - s + 1, T) + . . . + u n (l , T). (3.5) 

Each term u n (j, T), j = l — s,l — s + 1, . . . , /, on the right-hand side of relation (3.5) solves the following 
finite-difference counterpart to problem (2.14): 

^ ] ®mnttn(jt T) = fm(ji T), j ~ i S, / S -f- 1, . . . , /, 

" € ' V '" (3.6) 

w„(j, T) = 0, if n 4 r < ( j - 1 )T + r. 

The right-hand side f m (j, T) of (3.6) is given by the expression 


/ra(i, T) = fj(x,t,T) 
= 9j(t,T)f(x,t) 

(see (2.11), (2.12), (2.13)). Obviously, 


t—m4T 

x—(ni \ h, mzhMizh), t=m4T 


frn(j , T) = 0, if m 4 r < (j - 1)T. (3.7) 

Actually, it turns out that to obtain a discrete approximation (see (3.5)) to the solution Uj(x,t) of problem 
(2.14) for x € S(t), (l - 1)T < t < IT, one needs to solve the finite-difference problem (3.6) for every j, 
j = / - Sj l - s -f 1, . . . only on the time interval ( j - 1 )T < t < IT (rather than 0 < t < IT). The length 
of the largest of these intervals does not exceed d/(c - k) + 4 cT/(c - k) (see (3.3)) and does not depend on 
L Therefore, the following theorem holds. 

Theorem 3.1. The eiTor due to the replacement of the true solution Uj(x,t) in the exact formula (3.4) 
by the difference solution u\P = ^u n {j , T) of the first algorithm given by the approximate formula (3.5) 
does not exceed (s + 1 )Ch 2 , 


max 

ni ,n2,n 3 ,n 4 


u(nih, n>h, n 3 h, n 4 r) - u { J ] < (s + l)Ch 2 , 


(3.8) 


where the constant C depends only on the properties of the function f(x,t) for (l-s-l)T < t < IT and s+\ 
is the number of terms in the sum (3.5). Theorem 3.1 is, in fact, an obvious consequence of the construction 
of the first algorithm. Indeed, the estimate of type (1.12) will hold for each term on the right-hand side of 
formula (3.5). None of these terms is calculated on the time interval that exceeds T = d/(c-k) + 4 cT/(c - k) 
in length and therefore, all constants will be bounded. Consequently, the overall error estimate can be written 
in the form (3.8). Note, unlike in estimate (1.12), the constant C of (3.8) includes the factor of type \f\h.T 
that reflects the smoothness properties of the right-hand side. The function ^(t) of (2.9) that helps us build 
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the partition of unity (2.10) was chosen sufficiently smooth for the particular reason of guaranteeing that 
this partition that was introduced primarily for the algorithmic purposes will not “interfere” with the error 
estimates for the algorithm. Obviously, the error given by (3.8) will not increase, i.e., will remain uniformly 
bounded, when the number 714 of the time step increases as long as the smoothness properties of f(x^t) 
remain uniformly “good” with respect to t. Note, for higher-order schemes one may need smoother functions 
#(*)■ 

Let us now estimate the computer resources required by the proposed algorithm. Clearly, both the 
operations count and the amount of memory (i.e., number of words) needed for advancing one time step 
when calculating each term u n (j, T) of (3.5) are of the same asymptotic order 0(h ~ 3 ) with respect to the 
grid size h. The number of terms s 4- 1 does not change (i.e., does not grow) when the grid is refined as 
long as T is fixed. Therefore, neither does the overall number of arithmetic operations and words of memory 
required when calculating the solution by means of formula (3.5) - both quantities remain of the order 
0{h~ 3 ). The number of grid nodes that belong to the sphere S(t) for a fixed t - n 4 r is also of the order 
0 (h 3 ); therefore, the foregoing algorithm appears asymptotically non- improvable — linear with respect to 
the grid dimension — as long as one uses the scheme (1.6). 

We also note that out of the several terms that need to be computed according to (3.5), the first one, 
u n {l - 8 , T), is the most expensive numerically, its calculation by an explicit scheme up to the time level 
t — n 4 T — IT requires the widest grid domain of the size approximately 

/ d 4 cT \ 1 

d + 2(. S +l)Tc=d + 2(— + — ] - jc=—-^((3c-k)d + 8 c 2 T) (3.9) 

Dividing (3.9) by h and taking the third power of the result, we will obtain a quantity of the order 0(h~ 3 ) 
(as long as T is fixed while the grid is refined). This quantity gives an estimate of what will be the actual 
amount of memory needed for advancing one time step using the first algorithm. 

Finally, let us mention that when discussing the long-time integration we can consider a formulation that 
slightly differs from (3.1). Considering t from the interval (3.1) means, in fact, that l can be arbitrarily large 
and we calculate the solution on the interval of the fixed length T, which can be placed arbitrarily far away 
from the initial data. Alternatively, one may be interested in knowing the overall temporal evolution of the 
solution, i.e., in calculating the solution on an arbitrarily long time interval, say from 0 to some large Tfi na |. 
From the standpoint of building a lacunae-based algorithm that provides for non-accumulation of error and 
non-increasing expenses, this formulation is not much different from the one analyzed previously. For every 
time interval (3.1), T Tfi na i» he., every /, the solution will still be computed using formula (3.5). To 
advance further in time, we then need to replace l by l + 1 in formula (3.1). This will simply imply dropping 
the first term u n (l — s, T) on the right-hand side of formula (3.5) and adding the new last term u n (l 4- 1, T). 
In so doing, each term u(l — q , T), q — 0, 1, . . . , s — 1, for the previous interval l becomes u(l 4- 1 — (q 4- 1), T) 
for the new interval / 4 hi. Of course, for the new interval there is no need to calculate this term from the 
very beginning by solving the corresponding problem (3.6) starting from j = l - q; the calculation of each 
term of (3.5) that is not dropped when going from / to / 4- 1 (i.e., every term except the first one) is rather 
continued from the previous interval using the same explicit scheme. 

The Second Algorithm . In this algorithm, instead of formula (3.2) we use the following representation 
of the solution u{x,t) for x € S{t), (l - 1 )T < t < IT : 
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u{x,t) = vi- s (x,t,T,z) + vt-s+x(x,t,T,z) + . . . + vi(x,t,T,z), 


(3.10) 


here Vj(x, t, T, z), j = l — s,l — s 4- 1, . . . , /, are solutions of the corresponding problems (2.20). In so doing, 
the period 2 has to he chosen so that for every function uj(x , T), j = l - s, l - s + 1, . . . , /, the equality 


uj(x,t y T) = vj(x,t y T,z) 


hold on the entire time interval 


(j- DT<t< (j-l)T + 


d 2T (c -j- k ) 
c — k 


(3.11) 


on which according to formulae (2.17) the function Uj(x,t.,T) may differ from zero for x G S(t). In other 
words, we require that the time interval (3.11) be contained inside the time interval (2.22) (see Theorem 2.1) 
or 


W-i)r + 


d + 2T(c 4* k) 
• c — k 


<{j-l)T + 


z — d 
c + k y 


which yields the following condition: 


z > - 2 — (2c,d + 2(c + k) 2 T) . (3.12) 

C — rC 

Condition (3.12) essentially guarantees that all the waves generated by the sources x € S(t) that operate on 
the time interval 2 T will leave the domain of interest S(t) before the waves generated by the other sources 
from the periodic structure can enter this domain. 

To actually build the second algorithm, we replace the differential equation and initial condition of 
(2.20a) by the finite-difference equation and discrete initial condition, respectively: 


E UmnVnU, T, Z ) = f m (j, T, z), j = l - S, l - S + 1, . . . , l 

n£N m 

V n (j , T, z) — 0, if n 4 T < ( j - 1 )T + r, 

where the right-hand side f m {j, T, z) is a z-periodic grid function with node values 


(3.13a) 


fm(j, T, z) = fj(x,t,T,z)\ 

\x=(mih y m2h,,m-sh), f = m 4 r 

The periodic boundary conditions (2.20b) are replaced by their discrete counterparts in every coordinate 
direction jq , x>, and x^ (the period is 2 — rh): 
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V n (j, T, Z) = V n '(j, T, 2) 

n — (ni,n2,n 3 ,n 4 ), n = (ni -f Aqr, 712 4- 4- for, 7*4) (3.13b) 

fo , fo, fo — 0, ±1, ±2, . . . 

The approximation to the solution u(x,t) for x £ S(t ), (/ - 1)T < t < IT, in the second algorithm is 
obtained as an approximation to the sum (3.10) (compare to (3.5) and (3.2), respectively): 


u(ri\h , n^h, 773 / 7 , n 4 r) « u\l^ = 

T, z) + t7 n (/ — .s 4- 1, T, 2 ) + ... + !>«(/, T, 2 :). 


(3.14) 


Each term v n (j , T, z), j = / — 5 , / — 5 + 1, . . . , l , on the right-hand side of relation (3.14) solves the corre- 
sponding problem (3.13). Clearly, the error estimate of type (3.8) (with u\P of (3.5) replaced by uh 11 - of 
(3.14)) provided by Theorem 3.1 for the first algorithm, will hold for the second algorithm as well. 

As we did previously for the first algorithm, let us now consider the transition from / to / 4- 1 in formula 
(3.1) in the framework of the second algorithm. Assume we are interested in calculating the overall temporal 
evolution of the solution from t — 0 till an arbitrarily large t — 7k na ]. Over this period of time, the domain 
S(t) that was centered at x°(0) = rf(0), xf](0)) at the moment t = 0 can travel arbitrarily far in space 

from its initial location, in fact, as far as &7fi, ia i: 


ft dr 0 

x°(t) = (z°(t),x$(t),x$(t)) = x° (0) + J — -dt , (3.15a) 

\x°{t)\ < fcTfinai , 0 < t < T final . (3.15b) 

In the z-periodic setting, all functions are defined for \x a < zj 2, i — 1,2,3, and the edges jt, = ±z/ 2, 
i — 1,2,3, are identified with one another. Accordingly, instead of the motion described bv relation (3.15a), 
we consider the motion of S(t) on a three-dimensional toroidal surface. Instead of (3.15a) we will then have 

x°(t) = (x^(t),x? 2 (t),x%(t)) = l^(x°(0) + (3.16a) 

where {*} denotes the fractional part. Also, conforming to the periodicity conditions, inequality (3.15b) 
transforms into 


|*°(*)| < */2, 0 < t < Tfinal- (3.16b) 

The computational procedure does not change much. We calculate separately every term on the right-hand 
side of (3.14). When we go from l to / 4- 1, we stop calculating v n (l - s,T,z) and add the new term 
v n (l 4- I,T, z). This allows us to run the computation arbitrarily long with no error accumulation and 
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no growth of computer expenses per time step. In so doing, of course, the center x°(f) of 5(f), as well 
as the entire domain of interest 5(f) itself, can be located anywhere within the period, i.e., in the cube 
{l x i| 5- z l 2, i = 1,2,3}, and not necessarily close to its middle. This, however, does not affect the solution 
calculated inside 5(f) because according to condition (3.12), no waves from outside ran enter the domain 
S(t) before the waves generated inside 5(f) by the sources that operate during the interval 2 T leave it. Then, 
as soon as these waves have left, this entire portion of the solution, both the waves generated inside S(t ) 
that have already left it and the waves generated outside 5(f) by the other sources of the periodic structure 
that, operate during the same time interval it taken out by eliminating the term v n (l - s,T,z) in the sum 
(3.14). As has been pointed out, this does not essentially change anything inside S(t.) as the waves have 
already left, but it prevents the waves generated outside from propagating further in. 

The second algorithm appears somewhat more efficient than the first one as the actual domain size 
and consequently, the number of grid nodes involved, are smaller for the second algorithm, compare (3.12) 
against (3.9). The number of terms s + 1 in formulae (3.5) and (3.10) does not depend on the grid size h if 
T does not depend on h. If, however, T decreases as h decreases (when the grid is refined) then the number 
s and also constants involved in the estimates 1 0(h "*)| < K.h~ 3 grow. If T ~ t, then the number s becomes 
•s = 0(h~ l ) and the operations count, accordingly, 0(h~ 4 ) instead of the non-improvable quantity 0(h~ 3 ). 

The third algorithm keeps the non-improvable asymptotic order 0(h~ 3 ) for both the number of arith- 
metic operations and amount of memory required for advancing one time step even if T decreases when the 
grid is refined as long as T > r| In r|. 

The Third Algorithm. The third algorithm uses exactly the same approximation on each time step as 
the second algorithm does (see (3.10)). However, the computations in the third algorithm are organized in a 
substantially different way. In this algorithm, instead of calculating separately each term on the right-hand 
side of (3.14) we rather use a “one-sweep” time-marching approach and when it comes to the transition 
from / to / + 1 in (3.1), the term v„{l - s,T,z) on the right-hand side of (3.14) is taken out by the explicit 
subtraction. 

Introduce a new grid function 


^n(liSjT,z) — ^ ’ VnijiT.z), (3.17) 

j=l-s 

where v n (j,T,z) is the solution of problem (3.13). Here .s is the same integer number as in formula (3.10). 
Notice that for any n = (u ] h. n 2 h y n 3 h . n^T ) only several terms of the series (3.17) may differ from zero. For 
those nodes » that belong to the grid time levels 


t = ti^t, n-4 


(l-l)T (l-l)T IT 

_ 1 ' y • * • 1 


(3.18) 


the sum (3.17) coincides with the approximation of the solution u(x,t) given by formula (3.14) for 
(nih,n 2 h,n 3 h) e 5(f), t = n 4 r. Thus, the computation of V n {l,s,T,z) by formula (3.17) can be inter- 
preted as the computation of the approximate solution by formula (3.14). 

Substituting expression (3.17) into the left-hand side of the finite-difference equation (1.6), we obtain 
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(3.19) 


y: a mti v n (l, s, T, z) = /,„(/ - s, T, z) + 

n€N m 

+/m (t “ .S + l ,T,2) + ... + 2)-h ... 

According to the definition of $(£). see (2.9), and formula (2.11), the following equality holds 


/hi — / (rn 1 /i, TTi ‘2 h , m3 /1 , ^4 r ) — £ /m 0 \r,*) (3.20) 

*=*-• 

for all those m, for which 


m 4 r > (/ - 6*)T. (3.21) 

Equalities (3.19) and (3.20) imply that the function V n (l, 6, T, 2 ) satisfies the finite-difference equation 


r. a mn V„(l,s,T,z) = / m (3.22) 

n€A'm 

for those m, for which the inequality (3.21) holds. Clearly, for / = 0 and n 4 = 0, n 4 = 1, we have 
V n (l,s,T,z) = 0. 

Assume now that for some integer l > 0 we already know the values of V n (l, 6,T, z) on the first two time 
levels of the grid (3.18), i.e., for n 4 = (/ — l)T/r and n 4 = (/ — 1)T /r + 1. We will describe the computation of 
V n ( / , s , T, z) on all other levels of the grid (3.18), as well as computation of V n (l 4- 1, s, T, z) for n 4 = IT jr and 
TI 4 = IT /t 4- 1. Thus, we will have completed the cycle of going from the use of V n (l, 6, T, z) for calculating the 
approximate solution on the grid levels (3.18) to the use of V n (l 4- 1,6, T, z) for calculating the approximate 
solution on the subsequent levels 


t 


IT IT 

n 4 r, n 4 — — , hi, 

r r 


(/ 4- 1)T 
r 


- 1. 


So, having specified \ r n (l,s,T, z) for 


t = 7 } 4 r, n 4 


(*-pr 

T 


and 


(l - 1)T 

T 


+ It 


we calculate the values of V n (l, 6,T, z) on all other levels of the grid (3.18) using the explicit finite-difference 
scheme (3.22). Then, the values of V n (l 4- 1,6, T, z) for 714 = IT jr and n 4 = IT jr 4- 1 can be obtained with 
the help of the following obvious recurrence formula: 


V n (l 4- 1 ,s,T,z) - Vn(l, s,T,z) - v n (l - s,T,z). 


(3.23) 
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The first term on the right-hand side of formula (3.23) is calculated with the explicit scheme (3.22) using the 
data \ s,T, z) that have already been obtained on the last two levels 7*4 = IT jr — 2 and n\ ~ IT/r — 1 
of the grid (3.18). The second term v n (l - s , T, z) on the right-hand side of formula (3.23) is the solution of 
problem (3.13) for j = l — *#. The right-hand side f rn (j, T, z) of the latter problem may differ from zero only 
for those m = {mih.iibyh.rivji.m^r) that satisfy 


(l — S — 1 )T < 77l 4 T < (l — s - f 1)T, 

and the initial data for calculating v n (l — s,T,z) are homogeneous. 

The actual computation of v n (l — s,T , z) will be split into two stages. First, we will calculate the solution 
v n (l — s,T,z) with the explicit finite-difference scheme step after step in time for the levels 


t = 7*4 T, TI 4 


(l ~ « ~ 1 )T 


T 

The computation of v n (l - s, T, z) on the levels (3.24) takes 


(l-s + l)T 

T 


(3.24) 


2 T 

"T (ft) (3 - 25) 

arithmetic operations, where ft is the number of operations required for calculating the solution in one grid 
node on the next time level while the solution on the two previous levels is already known. (In other words, 
H is the number of arithmetic operations “on the stencil” of the scheme.) 

The values of the solution v n (l — s, T, z) on the last two levels of (3.24), i.e., n 4 = (l-s + 1 )T/r - 1 
and 7*4 = (/ — s + l)T/r, will be used as the initial data for calculating this solution v n (l - 5, T\ z) for 
7*4 = IT jr and 7*4 = IT /t + 1 . As the right-hand side f m (l - s , T, z) is equal to zero, f m (l - s, T, z) = 0, for 
77*4 > (/-.s-fl)T/r, an efficient w r ay to calculate the solution v n (l - a*,T, z) forn 4 = IT/r and n 4 =IT/t+ 1 
will be through representing it in the form of a discrete Fourier series while the initial data for v n (l — s,T,z) 
on the levels 7*4 = (/ — a + l)T/r — 1 and 7*4 = (/ — .s + l)T/r are known. The aforementioned finite Fourier 
expansion is built with respect to the following system of grid basis functions 


e J n = exp (i(n,j)), i = >/-!, 

n = (ni,n 2 ,n 3 ), j = (ji,j 2 ,h), 

ji,h, h = 0 , 1 , 2 , ... , \ - 1 , 
h 

(nj) = n x ji -j- n 2 h -f TI3J3- 


Hereafter we assume that zjh is a power of 2 so that the fast Fourier transform (FFT) can be used for for 
calculating the coefficients of the discrete Fourier series of a given grid function, as well as for restoring the 
point-wise values of the grid function from its Fourier representation. Each of the foregoing transformations 
requires 
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°(©M) 


(3.26) 


arithmetic operations. 

Thus, the second stage of the computation of v n (l - s,T, z) for n 4 = IT/t and n 4 = IT /r 4* 1 will first 
consist of Fourier transforming the data on the last two levels of grid (3.24), which takes (3.26) operations 
(twice the actual amount of (3.26) for two levels). Then, we advance each Fourier component independently 
to the levels n 4 = IT jr and n 4 = IT/r + 1 using the explicit formulae that are easily obtained from the 
Fourier representation of the finite-difference operator of (3.22); essentially, this reduces to multiplication by 
the appropriate powers of the corresponding amplification factors and obviously takes O ) 3 ^ operations. 
Finally, we need to restore the point-wise values of v n (l — s,T,z) at n 4 = IT/t and 77.4 = IT jr + 1 using the 
inverse FFT, which again takes (3.26) operations. The overall computational cost of this second stage will 
then be 


O (© 3| “0 + o ((i)V o ((I) 3l %)= c, (© , '"l) < 327 > 

arithmetic operations. Consequently, the total operations count for calculating v n (l — s, T, z) for 774 = IT/t 
and n 4 — IT/t- hi, i.e., calculating the second term in the recurrence formula (3.23), consists of the expenses 
for the first (3.25) and second (3.27) stages and adds up to 

< 3 - 28 > 

operations. Recurrence formula (3.23) is used for the entire “chunk” of 2 T/r time levels, (/ — s — 1)T < t < 
(l — s + 1)T. Therefore, if one recalculates the associated expense (3.28) proportionally per time step, it 
obviously reduces to 


The calculation of the first term in the recurrence formula (3.23) also requires O |) 3 ^) arithmetic operations 
per time step as this is done simply using the explicit scheme (3.22). 

Summarizing, we conclude that the overall computational cost of the third algorithm is O arith- 

metic operations per one time step. It is also easy to see that the required amount of memory (i.e., number 
of words) is of the order O as well. 

As the third algorithm essentially reproduces the calculation according to formula (3.14) with the differ- 
ence only in the computational procedure, the uniform error estimate of type (3.8) provided by Theorem 3.1 
for the first algorithm, will hold for the third algorithm as it does for the second algorithm. Moreover, 
the interpretation of spatial periodicity given by formulae (3.16) and subsequent comments for the second 
algorithm, applies to the third algorithm with no changes. 
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4. Possible Generalizations. First, let us note that the assumption of homogeneity of the initial data 
(1.2) can be alleviated and replaced by a weaker requirement 

= 

/=o 

where ipo{x) and <£>i(x) are sufficiently smooth functions that turn into zero outside the domain 5(t)| _ = 
5(0). 

Further, the requirement of smoothness for f(x,t) throughout the entire space-time (x, t) along with the 
consideration of /(x, t) only for t > 0 actually implies that /(x, t) and its first several derivatives with respect 
to t have to be smooth as t — > 4-0. This condition can also be alleviated by requiring that the function 
/(x,t), supp f{x,t) C S(t), be smooth for t > 0 rather than on the entire space-time (x, t). The resulting 
Cauchy problem, which appears somewhat more complex, can actually be reduced to problem (1.1), (1.2) if 
one represents the solution to the new problem as a sum of tw T o functions: 


u(x,t) 


t = o 


= ¥>o(s), 


du(x, t) 


dt 


u(x, t) = v(x, t) + w(x, t). 

The function u(x,t) will be a solution to the Cauchy problem with the given non- homogeneous initial data 
and the right-hand side F(x,t) = 9(t)f(x,t) that turns into zero for t > 1. The function w{x,t) will be a 
solution to the problem 


d 2 w 

W 


f d 2 w d 2 w d 2 w \ 

V dx'i + dx\ + dxj ) 


= f{x,t) -F(x,t), 


t > 0, 


W 


t - 0 


dw 

dt 


= 0 . 


(4.1) 


Problem (4.1) is obviously of the type (1.1), (1.2). The problem for v{x,t) needs to be solved only till some 
t — to, after which v(x,t) = 0 when x € S{t) because of the presence of lacunae in the solutions of the 
three-dimensional wave equation. 

Another possible generalization includes building similar lacunae-based algorithms for the long-time 
numerical integration of problems in acoustics (linearized Euler's equations), electromagnetics (Maxwell's 
equations), and elastodynamics (Lame’s equations). These algorithms may then be used for constructing 
highly accurate global genuinely time-dependent ABCs in the corresponding fields of scientific computing. 
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